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Abstract 

Most real-world networks are not isolated. In order to function fully, they are interconnected with other 
networks, and this interconnection influences their dynamic processes. For example, when the spread of a 
disease involves two species, the dynamics of the spread within each species (the contact network) differs 
from that of the spread between the two species (the interconnected network). We model two generic 
interconnected networks using two adjacency matrices, A and B, in which A is a 2N x 2N matrix that 
depicts the connectivity within each of two networks of size N, and B a 2N x 2N matrix that depicts 
the interconnections between the two. Using an N-intertwined mean-field approximation, we determine 
that a critical susceptable-infected-susceptable (SIS) epidemic threshold in two interconnected networks 
is l/\i(A + aB), where the infection rate is /3 within each of the two individual networks and a/3 in 
the interconnected links between the two networks and Xi(A + aB) is the largest eigenvalue of the matrix 
A + aB. In order to determine how the epidemic threshold is dependent upon the structure of interconnected 
networks, we analytically derive Xi(A + aB) using perturbation approximation for small and large a, the 
lower and upper bound for any a as a function of the adjacency matrix of the two individual networks, and the 
interconnections between the two and their largest eigenvalues/eigenvectors. We verify these approximation 
and boundary values for Xi(A + aB) using numerical simulations, and determine how component network 
features affect X\(A + aB). We note that, given two isolated networks Gi and G2 with principle eigenvectors 
x and y respectively, Ai (A + aB) tends to be higher when nodes i and j with a higher eigenvector component 
product Xiyj are interconnected. This finding suggests essential insights into ways of designing interconnected 
networks to be robust against epidemics. 



1 Introduction 

Complex network studies have traditionally focused on single networks in which nodes represent agents and links 
represent the connections between agents. Recent efforts have focused on complex systems that are comprised 
of interconnected networks, a configuration that more accurately represents real- world networks [HE]. Real- 
world power grids, for example, are almost always coupled with communication networks. Power stations need 
communication nodes for control and communication nodes need power stations for electricity. The influence 
of coupled networks on cascading failures has been widely studied [IJ 02 [H [5l |6] . When a node at one end of an 
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interdependent link fails, the node at the other end of the link usually fails. A non-consensus opinion model of 
two interconnected networks that allows the opinion interaction rules within each individual network to differ 
from those between the networks was recently studied J7J. This model shows that opinion interactions between 
networks can transform non-consensus opinion behavior into consensus opinion behavior. 

In this paper we investigate the susceptable-infected-susceptable (SIS) behavior of a spreading virus, a 
dynamic process in interconnected networks that has received significant recent attention [8j [9j [101 [H] . An 
interconnected networks scenario is essential when modeling epidemics because diseases spread across multiple 
networks, e.g., across multiple species or communities, through both contact network links within each species 
or community and interconnected network links between them. Dickison et al. 9 study the behavior of 
susceptible-infected-recovered (SIR) epidemics in interconnected networks. Depending on the infection rate in 
weakly and strongly coupled network systems, where each individual network follows the configuration model 
and interconnections are randomly placed, epidemics will infect none, one, or both networks of a two-network 
system. Mendiola et al. [10] show that in SIS model an endemic state may appear in the coupled networks even 
when an epidemic is unable to propagate in each network separately. In this work we will explore how both 
the structural properties of each individual network and the behavior of the interconnections between them 
determine the epidemic threshold of two generic interconnected networks. 

In order to represent two generic interconnected networks, we represent a network G with N nodes using an 
N x N adjacency matrix A\ that consists of elements ay, which are either one or zero depending on whether 
there is a link between nodes i and j. For the interconnected networks, we consider two individual networks G\ 
and G2 of the same size N. When nodes in G\ are labeled from 1 to N and in G2 labeled from N + 1 to 2N, 
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the two isolated networks G\ and G2 can be presented by a 2N x 2N matrix A = 



A 2 

their corresponding adjacency matrix A\ and A2 respectively. Similarly, a 2N x 2N matrix B = 



composed of 



B 12 

bL 

represents the symmetric interconnections between G% and G2. The interconnected networks are composed of 
three network components: network A\, network A2, and interconnecting network B. 

In the SIS model, the state of each agent at time t is a Bernoulli random variable, where Xi(t) = if node i is 
susceptible and Xi(t) — 1 if it is infected. The recovery (curing) process of each infected node is an independent 
Poisson process with a recovery rate 6. Each infected agent infects each of its susceptible neighbors with a rate 
j3, which is also an independent Poisson process. The ratio r = j3/S is the effective infection rate. A phase 
transition has been observed around a critical point t c in a single network. When t > r c , a non-zero fraction of 
agents will be infected in the steady state, whereas if t < t c , infection rapidly disappears |12|I13]. The epidemic 
threshold via the N-intertwined mean- field approximation (NIMFA) is t c = Xl ^ A ^ , where X±(A) is the largest 
eigenvalue of the adjacency matrix, also called the spectral radius [14]. For interconnected networks, we assume 
that the curing rate 6 is the same for all the nodes, that the infection rate along each link of G\ and G2 is (3, 
and that the infection rate along each interconnecting link between G\ and G2 is a/3, where a is a real constant 
ranging within [0, 00) without losing generality. 

We first show that the epidemic threshold for (3/5 in interconnected networks via NIMFA is r c = ^TTXhxB) ' 
where \\(A + aB) is the largest eigenvalue of the matrix A + aB. We further express \i(A + aB) as a function 
of each network component A\ 1 A2, and B and their eigenvalues/eigenvectors to reveal the contribution of 
each component network. This is a significant mathematical challenge, except for special cases, e.g., when 
A and B commute, i.e., AB = BA (see Sec. 13. ip . Our main contribution is that we analytically derive for 
the epidemic characterizer Xi(A + aB) (a) its perturbation approximation for small a, (b) its perturbation 
approximation for large a, and (c) its lower and upper bound for any a as a function of component network 
Ai, A 2 , and B and their the largest eigenvalues/eigenvectors. Numerical simulations in Sec. 0] verify that these 
approximations and bounds well approximate Xi(A + aB), and thus reveal the effect of component network 
features on the epidemic threshold of the whole system of interconnected networks, which provides essential 



2 



insights into designing interconnected networks that are robust against the spread of epidemics (see Sec. [SJ. 

Sahnch ct al. TT] recently studied SIS epidemics on generic interconnected networks in which the infection 
rate can differ between G\ and G2, and derived the epidemic threshold for the infection rate in one network 
while assuming that the infection does not survive in the other. Their epidemic threshold was expressed as 
the largest eigenvalue of a function of matrices. Our work explains how the epidemic threshold of generic 
interconnected networks is related to the properties (eigenvalue/eigenvector) of each network component A\, 
A2, and B without any approximation on the network topology. 

Graph spectra theory |15j and modern network theory, integrated with dynamic systems theory, can be used 
to understand how network topology can predict these dynamic processes. Youssef and Scoglio [16] have shown 
that a SIR epidemic threshold via NIMFA also equals 1/Ai. The Kuramoto synchronization process of coupled 
oscillators [17] and percolation [18] also features a phase transition that specifies the onset of a remaining 
fraction of locked oscillators and the appearance of a giant component, respectively. Note that a mean-field 
approximation predicts both phase transitions at a critical point that is proportional to 1/Aj.. Thus we expect 
our results to apply to a wider range of dynamic processes in interconnected networks. 



2 Epidemic Threshold of Interconnected Networks 

In the SIS epidemic spreading process, the probability of infection = E[Xi(t)] for a node i in interconnected 
networks G is described by 



dvjjt) 
dt 



2N 2N 

p J2 a ijVj (i) + ap Yj bi M W ] 0- - v o (*)) ~ 5v J (*) 
j=i 3=1 



via NIMFA, where ay and 6^ is an element of matrix A and B respectively. Its matrix form becomes 

= (/3(A + aB) V(t) ~ SI) - Pdiag ( Wj (t)) (A + aB) V{t). 
The governing equation of the SIS spreading process on a single network A\ is 

= (pAtV® - SI) - Pdiag (v t (t)) A ± V(t), 

at 

whose epidemic threshold has been proven [14] to be 

1 



Ai(Ax)' 



which is a lower bound of the epidemic threshold [19] . Hence, the epidemic threshold of interconnected networks 
by NIMFA is 

(i) 



Xi(A + aB) 

which depends on the largest eigenvalue of the matrix A + aB. The matrix A + aB is a weighted matrix, where 
< a < 00. Note that the NIMFA model is an improvement over earlier epidemic models [13] in that it applies 
no approximations to network topologies, and thus it allows us to identify the specific role of a general network 
structure on the spreading process. 

3 Analytic approach: \i(A + aB) in relation to component network 
properties 

The spectral radius Xi(A + aB) as shown in the last section is able to characterize epidemic spreading in 
interconnected networks. In this section we explore how Xi(A + aB) is influenced by the structural properties 



3 



of interconnected networks and by the relative infection rate a along the interconnection links. Specifically, we 
express \i(A + aB) as a function of the component network A 1: A 2 , and B and their eigenvalues/eigenvectors. 
(For proofs of theorems or lemma, see the Appendix.) 



3.1 Special cases 



We start with some basic properties related to Ai (A + aB) and examine several special cases in which the 
relation between Xx(A + aB) and the structural properties of network components A\, A 2 and B are analytically 
tractable. 

The spectral radius of a sub-network is always smaller or equal to that of the whole network. Hence, 



Lemma 1 



Lemma 2 



Xx{A + aB) > Ai(A) =max(A 1 (A 1 ), Xx(A 2 )) 



Ai (A + aB) > aXi(B) 
The interconnection network B forms a bipartite graph. 



Lemma 3 The largest eigenvalue of a bipartite graph B = 
where Bxi is possibly asymmetric \1 5]j . 



B 12 

2 



Bio 



follows Ai(B) = JX X (B( 2 B 12 



Lemma 4 When G\ and G 2 are both regular graphs with the same average degree E[D] and when any two 
nodes from G\ and G 2 respectively are randomly interconnected with probability pi, the average spectral radius 
of the interconnected networks follows 

E[Xx(A + aB)] = E[D] + aN Pl 
if the interdependent connections are not sparse. 

A dense Erdos-Renyi (ER) random network approaches a regular network when N is large. Lemma 01 thus, 
can be applied as well to cases where both G\ and G 2 are dense ER random networks. 

For any two commuting matrices A and B, thus AB = BA, X±(A + B) — X±(A) + X±(B) [15]. This property 
of commuting matrices makes the following two special cases analytically tractable. 



Lemma 5 When A + aB — 



aB = 


' Ax 
















+ a 











At _ 










, i.e., the interconnected networks are composed of 

two identical networks, where one network is indexed from 1 to N and the other from N + 1 to 2N , with 
an interconnecting link between each so-called image node pair (i, N + i) from the two individual networks 
respectively, its largest eigenvalue Xi(A + aB) = X\(A) + a. 



Proof. When A + aB 



' A, 





+ a 





/ " 





Ax _ 




/ 






matrix A and aB are commuting 



A - aB = a 



Ax 
Ax 



= aBA 



Therefore, Xx{A + aB) — Xx(A) + Xx(aB) = Xi(Ai) + aXi(B). The network B is actually a set of isolated links. 
Hence, Xx(B) = 1. ■ 



Lemma 6 When A + aB — 



' Ax 





+ a 





Ax ' 





Ax _ 










, its largest eigenvalue Xx(A + aB) = (1 + a) Xx(Ax). 
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Proof. When A + aB — 



' A\ 










A x ' 




+ a 









A x _ 










matrix A and aB are commuting 



A - aB = a 



A\ 
A\ 



aBA 



Therefore Xx(A + aB) = Xx(A) + Xx(aB) = (1 + a) Xx(A) = (1 + a) Xx(A x ). ■ 

When A and B are not commuting, little can be known about the eigenvalues of Ai(yl + aB), given the 
spectrum of A and of B. For example, even when the eigenvalue of A and B are known and bounded, the 
largest eigenvalue of Ai(^4 + aB) can be unbounded [15] . 

3.2 Lower bounds for Ai (A + aB) 

We now denote matrix A + aB to be W. Applying the Rayleigh inequality [T51 p. 223] to the symmetric matrix 
W = A + aB yields 

z T Wz , ,„„ 
— ~— < Xx (W) 

z 1 z 

where equality holds only if z is the principal eigenvector of W. 

Theorem 7 The best possible lower bound z z ^ z z of interdependent networks W by choosing z as the linear 
combination of x and y, the largest eigenvector of A\ and A2 respectively, is 



Xx(W)>max(Xx (Ax),Xx (A 2 )) + 
where £ = ax T B\2y. 



Xx (Ax) - Xx (A 2 ) 



e- 



Xx (Ax) - Ai (A 2 ) 



(2) 



When a = 0, the lower bound becomes the exact solution Ai (W) — Xl- When the two individual networks 
have the same largest eigenvalue Ai (Ax) — Xx (A 2 ), we have 

Ai (W)>Xx (Ax) + ax T B 12 y 

Theorem 8 The best possible lower bound A^ (W) > z z by choosing z as the linear combination of x and 
y, the largest eigenvector of Ax and A 2 respectively, is 

(A? (Ax) + a 2 \\Bf 2 x\\ 2 2 + X\ (A 2 ) + a 2 \\B 12 y 



X( (W) > 



(3) 



\ 



'A 2 (Ax) + a 2 \\B{ 2 x\\l - X\ (A 2 ) - a 2 \\B 12 y\\ 



2 



where 9 = a(X x (Ax) + Xx (A 2 )) x T B 12 y. 
In general, 

The largest eigenvalue is lower bounded by 



z T W k z 



z T W k z 



< X k (W) 

1/k 

< Xx (W) 



/T s\l/ S / T k \ / ^ 

Theorem 9 Given a vector z, z ^ z < z ^ z when k is an even integer and < s < k. Further- 
more, 

' z T W k z" 



limk- 



z 1 z 



= Ai (W) 
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Hence, given a vector z, we could further improve the lower bound ( z Jr z z ) by taking a higher even 
power k. Note that Theorem [7] and [8] express the lower bound as a function of component network A\, A 2 and B 
and their eigenvalues/eigenvectors, which illustrates the effect of component network features on the epidemic 
characterizer Ai (W). 

3.3 Upper bound for A x (A + cuB) 

Theorem 10 The largest eigenvalue of interdependent networks \\{W) is upper bounded by 

Ai (W) < max (Ai {A x ),X 1 (A 2 )) + aX x (B) (4) 



= max(Ai(A 1 ),Ai(A 2 ))+a v /Ai(Bi 2 B 1 J 2 ) (5) 

This upper bound is reached when the principal eigenvector of B\ 2 B\ 2 coincides with the principal eigenvector 
of A\ if Ai(Ai) > A 1 (A 2 ) and when the principal eigenvector of Bf 2 Bi 2 coincides with the principal eigenvector 
of A 2 if Ai(Ax) < X 1 (A 2 ). 

3.4 Perturbation analysis for small and large a 

In this subsection, we derive the perturbation approximation of \i(W) for small and large a, respectively, as a 
function of component networks and their eigenvalues/eigenvectors. 

We start with small a cases. The problem is to find the largest eigenvalue sup z ^ ' "' Jr* °f W) with the 
condition that 

f (W-XI)z = 
\ z T z = l 

When the solution is analytical in a, we express A and z by Taylor expansion as 

oo 

A = ]TA( fe >a fe 

oo 
k=0 

Substituting the expansion in the eigenvalue equation gives 



(A + aB) z (k) a k = A(fc) « fc E z(k)ak 

k=0 k=0 k=0 

where all the coefficients of a k on the left must equal those on the right. Performing the products and reordering 
the series we obtain 



jr (AzW + B Z (*-D) a k = ]T j £ A( fe - 



k=0 k=Q \i=0 

This leads to a hierarchy of equations 

k 



k-i) z (i) 



i=0 

The same expansion must meet the normalization condition 



z T z = 1 



or equivalently, 



£z«>a* ] =1 
fc=o j=o 
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where (w, v) — ^ ■ u^Vi represents the scalar product. The normalization condition leads to a set of equations 



u 

i=l 







(0) 



for any k > 1 and (z* ', z' ') = 1. 

Let Xi(Ai) (Xi(A 2 )) and denote the largest eigenvalue and the corresponding eigenvector of Ai(A 2 ) 
respectively. We examine two possible cases: (a) the non-degenerate case when Xi(Ai) > Xi(A 2 ) and (b) the 
degenerate case when X\(Ai) = X\(A 2 ) and the case Xi(Ai) < X\(A 2 ) is equivalent to the first. 



Theorem 11 For small a, in the non-degenerate case, thus when Xi(A\) > Xi(A 2 ), 

Ai(W) = XxiAx) + a 2 {x^) T B 12 [X X {A X )I - A 2 )^ B&x® + 0(a 3 ) 
where (z^) T = x T T ) . 



(7) 



Note that in (fT4)l B is symmetric and (XS^I — A) is positive definite and so is B(\Wl-A) B. Hence, 
this second order correction A^ 2 -* is always positive. 

Theorem 12 For small a, when the two component networks have the same largest eigenvalue Xi(Ai) = 
Xi(A 2 ), 



X 1 {W) = X X {A X ) + \ax T B X2 y + c?{y^) T Bf 2 (X^I - A l+ ^(#) T )- 1 - 



'l2\ /y ' ± ~ ^1 ~r A ' ' \^ •■) ) ■ (8) 
(B 12 y^ - \(Q x <f>) + ( x (°yfB 12 (X^I -A 2 + x^ix^fr'B^ ~ \WyW) + 0(a 3 ) 



' A, 





and B = 





At ' 













At _ 










In the degenerate case, the first order correction is positive and the slope depends on B\ 2l y, and x. When 
Ai and A 2 are identical, the largest eigenvalue of the interdependent networks becomes 

A = Ai(Ai) + a (B 12 x,x) + 0{a 2 ) 

When A = and B = , our result (HI) in the degenerate case up to the first order leads to 

Ai 10 

X\(A+aB) = Xi(A)+a, which is an alternate proof of Lemma[5] When A = 
(J5J) again explains Lemma [5] that X%(A + aB) = (1 + a) Ai(Ai). 

Lemma 13 For large a, the spectral radius of interconnected networks is 

Xi(A + aB) = aXi(B) + v T 'Av + O (a" 1 ) (9) 
where v is the eigenvector belonging to X\(B) and 

Xi(A + aB) < Xi(A) + aXi(B) + O (a" 1 ) 

Proof. The lemma [13] follows by applying perturbation theory [20] to the matrix a(B + —A) and the Rayleigh 
principle j!5) . which states that v T Av < Ai(A), for any normalized vector v such that v T v — 1, with equality 
only if v is the eigenvector belonging to the eigenvalue X\{A). m 
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4 Numerical simulations 



In this section, we employ numerical calculations to quantify to what extent the perturbation approximation 
([7J and ([8|) for small a, the perturbation approximation (j9]) for large a, the upper (j4]) and lower bound ((3]) are 
close to the exact value Xi(W) = Xi(A + aB). We investigate the condition under which the approximations 
provide better estimates. The analytical results derived earlier are valid for arbitrary interconnected network 
structures. For simulations, we consider two classic network models as possible topologies of G\ and G2' (i) 
the Erdos-Renyi (ER) random network [21] [22l E2] and (ii) the Barabasi- Albert (BA) scale- free network [24] . 
ER networks are characterized by a Binomial degree distribution Pr[D = k] = ( N ~ 1 ^p k (l — p) JV ~ 1 ~ fc , where 
N is the size of the network and p is the probability that each node pair is randomly connected. In scale-free 
networks, the degree distribution is given by a power law Pr[Z) = k] = ck~ x such that X^fe!] 1 C ^~ A = 1 & n d 
A = 3 in BA scale-free networks. 

In numerical simulations, we consider N\ = N2 = 1000. Specifically, in the BA scale-free networks m = 3 
and the corresponding link density is pba — 0.006. We consider ER networks with the same link density 
Per = Pba — 0.006. A coupled network G is the union of G\ and G2, which are chosen from the above 
mentioned models, together with random interconnection links with density pi, the probability that any two 
nodes from G\ and G2 respectively are interconnected. Given the network models of G± and G2 and the 
interacting link density pi, we generate 100 interconnected network realizations. For each realization, we 
compute the spectral radius X\(W), its perturbation approximation ([7]) and (JHJ for small a, the perturbation 
approximation ([9]) for large a, upper bound ((4]) and lower bound <j3j> for any a. We compare their averages over 
the 100 coupled network realizations. We investigate the degenerate case Ai(Gi) = \\{G2) where the largest 
eigenvalue of G\ and G2 are the same and the non-degenerate case where Ai(Gi) ^ Ai(G2) respectively. 



4.1 Non-degenerate case 
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(a) *) 
Figure 1: A plot of \i(W) as a function of a for both simulation results (symbol) and its (a) perturbation 
approximation ([7} for small a (dashed line) and (b) perturbation approximation ([9]) for large a (dashed line). 
The interconneted network is composed of an ER random network and a BA scale-free network both with 
N = 1000 and link density p — 0.006, randomly interconnected with density pi. All the results are averages of 
100 realizations. 

We consider the non-degenerate case in which G\ is a BA scale-free network with N = 1000, m — 3, G2 is an 
ER random network with the same size and link density per — Pba — 0.006, and the two networks are randomly 
interconnected with link density pi. We compute the largest average eigenvalue i?[Ai(W / )] and the average of 
the perturbation approximations and bounds mentioned above over 100 interconnected network realizations for 
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each interconnection link density pi € [0.00025,0.004] such that the average number of interdependent links 
ranges from -j, ¥-,N, 2N to 4JV and for each value a that ranges from to 10 with step size 0.05. 

For a single BA scale- free network, where the power exponent (3 — 3 > 2.5, the largest eigenvalue is 
(1 + o(l)) %/^max where d max is the maximum degree in the network [25]. The spectral radius of a single ER 
random graph is close to the average degree (N — 1)per. when the network is not sparse. When pi = 0, 
Ai(G) = max (Xi(Ger), Xx(Gba)) — X\(Gba) > Xi(Ger)- The perturbation approximation is expected to 
be close to the exact Xi(W) only for a — > and a — > oo. However, as shown in Fig. Ufa), the perturbation 
approximation for small a approximates \±(W) well for a relative large range of a, especially for sparser 
interconnections, i.e., for a smaller interconnection density p\. Figure [TJb) shows that the exact spectral radius 
\\{W) is already close to the large a perturbation approximation, at least for a > 8. 
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(a) 0») 
Figure 2: Plot X\(W) as a function of a for both simulation results (symbol) and its (a) its lower bound ((3]) 
(dashed line) and (b) upper bound ^ (dashed line) . The interconneted network is composed of an ER random 
network and a BA scale-free network both with TV = 1000 and link density p = 0.006, randomly interconnected 
with density pj. All the results are averages of 100 realizations. 

As depicted in Fig. [51 the lower bound © and upper bound (Q} are sharp, i.e., close to Ai (W) for small a. The 
lower and upper bounds are the same as Ai(W) when a — > 0. For large a, the lower bound better approximates 
Ai(VF) when the interconnections are sparser. Another lower bound aXi(B) < Ai(W), i.e., Lemma [2 is sharp 
for large a, as shown in Fig. [3J especially for sparse interconnections. We do not illustrate the lower bound ([2]) 
because the lower bound (j3]) is always sharper or equally good. The lower bound aAi(-B) considers only the 
largest eigenvalue of the interconnection network B and ignores the two individual networks G\ and Gi. The 
difference Ai(W)— a\i(B) = v T Av + O (a^ 1 ) according to the large a perturbation approximation, is shown 
in Fig. |3] to be larger for denser interconnections. It suggests that G\ and Gi contribute more to the spectral 
radius of the interconnected networks when the interconnections are denser in this non-degenerate case. For 
large a, the upper bound is sharper when the interconnections are denser or when pi is larger, as depicted in 
FigureJ^fb). This is because aXi(B) < Xi(W) < aXi(B) + max(Ai(Ax), Xx(A2)). When the interconnections 
are sparse, X\(W) is close to the lower bound aXi(B) and hence far from the upper bound. 

Most interdependent or coupled networks studied so far assume that both individual networks have the same 
number of nodes N and that the two networks are interconnected randomly by N one-to-one interconnections, or 
by a fraction q of the N one-to-one interconnections where 0<g<l[TJ[Sl[7]- These coupled networks correspond 
to our sparse interconnection cases where pi < 1, when Xi(B) is well approximated by the perturbation 
approximation for both small and large a. The spectral radius X\(W) increases quadratically with a for 
small a, as described by the small a perturbation approximation. The increase accelerates as a increases and 
converges to a linear increase with a, with slope Xi(B). Here we show the cases in which Gi, G2, and the 
interconnections are sparse, as in most real-world networks. However, all the analytical results can be applied 
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a 

Figure 3: Plot \i(W) as a function of a for both simulation results (symbol) and its lower bound aXi(B) 
(dashed line). The interconneted network is composed of an ER random network and a BA scale- free network 
both with N = 1000 and link density p — 0.006, randomly interconnected with density pj. All the results are 
averages of 100 realizations. 

to arbitrary interconnected network structures. 



4.2 Degenerate case 
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(b) 

Figure 4: A plot of Xi(W) as a function of a for both simulation results (symbol) and its (a) perturbation 
approximation ([8]) for small a (dashed line) and (b) perturbation approximation ([9]) for large a (dashed line). 
The interconneted network is composed of two identical BA scale-free networks with N = 1000 and link density 
p = 0.006, randomly interconnected with density pj. All the results are averages of 100 realizations. 

We assume the spectrum |26j to be an unique fingerprint of a large network. Two large networks of the 
same size seldom have the same largest eigenvalue. Hence, most interconnected networks belong to the non- 
degenerate case. Degenerate cases mostly occur when G\ and Gi are identical, or when they are both regular 
networks with the same degree. We consider two degenerate cases where both network G\ and G2 are ER 
random networks or BA scale-free networks. Both ER and BA networks lead to the same observations. Hence 
as an example we show the case in which both G\ and G2 are BA scale-free networks of size N = 1000 and 
both are randomly interconnected with density pi € [0.00025, 0.004], as in the non-degenerate case. Figure @|a) 
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shows that the perturbation analysis well approximates Xi(W) for small a, especially when the interconnection 
density is small. Moreover, the small a perturbation approximation performs better in the non-degenerate case, 
i.e., is closer to Ai(W) than in degenerate cases [see Fig. [Ha)]. Similarly, Fig. [5] shows that both the lower and 
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(a) 00 
Figure 5: Plot Xi(W) as a function of a for both simulation results (symbol) and its (a) its lower bound ([3]) 
(dashed line) and (b) upper bound (HI (dashed line). The interconneted network is composed of two identical 
BA scale-free networks N — 1000 and link density p = 0.006, randomly interconnected with density pi. All the 
results are averages of 100 realizations. 

upper bound are sharper for small a. The lower bound better approximates Ai(W) for sparser interconnections 
whereas the upper bound better approximates Xi(W) for denser interconnections. 

Thus far we have examined the cases where Gi, G2, and the interconnections are sparse, as is the case 
in most real-world networks. However, if both G\ and G2 are dense ER random networks and if the random 
interconnections are also dense, the upper bound is equal to X\(W), i.e., X\{W) = Ai(Gi) + aXi(B) (see 
Lemma HJ. Equivalently, the difference Ai(W) — aXi(B) is a constant Ai(Gi) = Ai(G2) independent of the 
interconnection density pi . 

In both the non-degenerate and degenerate case, Ai(MK) is well approximated by a perturbation analysis for 
a large range of small a, especially when the interconnections are sparse, and also for a large range of large a. 
The lower bound (J3j) and upper bound Q are sharper for small a. Most real-world interconnected networks 
are sparse and non-degenerate, where our perturbation approximations are precise for a large range of a, and 
thus reveal well the effect of component network structures on the epidemic characterizer Xx(W). 



5 Conclusion 

We study interconnected networks that are composed of two individual networks G\ and G2 , and interconnecting 
links represented by adjacency matricies A\, A 2 , and B respectively. We consider SIS epidemic spreading in 
these generic coupled networks, where the infection rate within G\ and G2 is /?, the infection rate between the two 
networks is af3, and the recovery rate is 5 for all agents. Using a NIMFA we show that the epidemic threshold 

Ai 
A 2 

networks G\ and G2. The largest eigenvalue X\{A + aB) can thus be used to characterize epidemic spreading. 
This eigenvalue Ai (A + aB) of a function of matrices seldom gives the contribution of each component network. 
We analytically express the perturbation approximation for small and large a, lower and upper bounds for any 
a, of Ai(^4 + Q!-B) as a function of component networks Ay, A2, and B and their largest eigenvalues/eigenvectors. 
Using numerical simulations, we verify that these approximations or bounds approximate well the exact Ai(^4 + 
aB), especially when the interconnections are sparse and when the largest eigenvalues of the two networks G\ 



with respect to j3/6 is r c = XxfA+aB) > w ^ ere A = 



is the adjacency matrix of the two isolated 
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and Gi arc different (the non-degenerate case), as is the case in most real- world interconnected networks. Hence, 
these approximations and bounds reveal how component network properties affect the epidemic characterizer 
\\{_A + aB). Note that the term x T Bi 2 y contributes positively to the perturbation approximation flSJ) and 
the lower bound ^ of Xi(A + aB) where x and y are the principal eigenvector of network G\ and G 2 . This 
suggests that, given two isolated networks G\ and G 2 , the interconnected networks have a larger \\{A + aB) 
or a smaller epidemic threshold if the two nodes i and j with a larger eigenvector component product Xiyj from 
the two networks, respectively, are interconnected. This observation provides essential insights useful when 
designing interconnected networks to be robust against epidemics. The largest eigenvalue also characterizes the 
phase transition of coupled oscillators and percolation. Our results apply to arbitrary interconnected network 
structures and are expected to apply to a wider range of dynamic processes. 
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A Proofs 

A.l Proof of Lemma [4] 

In any regular graph, the minimal and maximal node strength are both equal to the average node strength. 
Since the largest eigenvalue is lower bounded by the average node strength and upper bounded by the maximal 
node strength as proved below in Lemma ll4[ a regular graph has the minimal possible spectral radius, which 
equals the average node strength. When the interdependent links are randomly connected with link density 
pi, the coupled network is asymptotically a regular graph with average node strength E[D] + aNpi, if pi is a 
constant. 

Lemma 14 For any N x N weighted symmetric matrix W , 

E[S] < Xi(W) < maxs r 

where s r = X)^=i w rj * s defined as the node strength of node r and E[S] is the average node strength over all 
the nodes in graph G. 

Proof. The largest eigenvalue Ai follows 

. x T Wx 

\\ = SUp Tf, 

x^O X 1 X 

when matrix W is symmetric and the maximum is attained if and only if x is the eigenvector of W belonging 



to \i{W). For any other vector y 7^ x, it holds that Ai > v y T^ v ■ By choosing the vector y = u = (1, 1, 1), 
we have 



iV iV N 
i— 1 j — 1 i—1 

where lOy is the element in matrix W and E[S] is the average node strength of the graph G. The upper bound 
is proved by Gerschgorin circle theorem. Suppose component r of eigenvector x has the largest modulus. The 
eigenvector can be always normalized such that 

x\ x 2 x r -i x r+ i xn 
, , . . . , , 1, , . . . , 

X ' f X y OCj' Xjy 3jf 
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where 



< 1 for all j. Equating component r on both sides of the eigenvalue equation Wx = \\x gives 



N N 
0=1 r 3=1 



N 

< J2\w rj \ 



when none of the elements of matrix W are negative. Since any component of x may have the largest modulus, 
Ai(W) < maxs r . ■ 



A. 2 Proof of Theorem 

We consider the 2N x 1 vector z as z T 



the linear combination of the principal eigenvector 



Cix T C 2 y T 

x and y of the two individual networks respectively, where x T x — 1, y T y = 1, Cj* + Cf = 1 such that z T z = 1 
and compute 



C l2 ; T CV 



Ax «-Bi2 Cix 
aBf 2 A 2 [C 2 y 

= Cfx T A lX + C 2 y T A 2 y + 2a2C 1 C 2 x T B 12 y 

= C 2 1 X 1 (A x ) + Ai (A 2 ) + 2dC 2 £ 

where £ = ax T B\ 2 y. By Rayleigh's principle Ai (W) > z Jfj' — z T Wz. We could improve this lower bound 
by selecting z as the best linear combination (C\ and C 2 ) of x and y. Let Ai be the best possible lower bound 

z T Wz 



7T- 



via the optimal linear combination of x and y. Thus, 



X L = max C 2 1 \ 1 {A 1 ) + Cl\ 1 {A 2 ) + 2C 1 C 2 i 

cf+cl=i 

We use the Lagrange multipliers method and define the Lagrange function as 

A = C 2 Ai [Ax) + C 2 2 Ai (A 2 ) + 2CxC 2 £ - fi (C? + C\ - l) 
where \i is the Lagrange multiplier. The maximum is achieved at the solutions of 

= 2CiAi (Ax) + 2C 2 £ - 2C 1M = 
= 2C 2 Ai (A 2 ) + 2CiC - 2C 2Ai = 

ou 2 

^ = C 2 + C 2 -1 = 

Note that fCxJ^j- + C 2 ^4j^ /2 = Al — ^ = 0, which leads to /i = A^. Hence, the maximum Al is achieved at 
the solution of 

C x Ai (A^ + C^-dXL = 
C 2 Aj (A 2 ) + Cxt - C 2 X L = 



that is 



det 



Ai (A x ) - A, 



Ai (A 2 ) - A z 



= 



13 



This leads to 



v _ Ai (Ax) + Ai (A 2 ) | J^ \ 1 (A 1 )^\ 1 (A 2 y ] ;1 



\i (A,) + Xx (A 2 ) 



Ai (Ax) - Ai (A 2 ) 



max(Ai (Ai),Ai (A 2 )) + 




Ai (Ai) - Ax (A 2 ) 



'Ai (^i)-Ai (A 2 ) 



Ax (At) - Ai (A 2 ) 



Ai (Ax) ~ Ai (A 2 ) 



The maximum is obtained when 



z T = ± 



Xi(A 2 )-X L T 
X 1 (A 1 )+\ 1 (A 2 )-2X L Jj 



X 1 (A 1 )+\ 1 (A 2 )-2X L y 



A. 3 Proof of Theorem H 



By Rayleigh's principle X\(W)> 
of x and y. The lower bound 



z T W 2 z. We consider z as linear combination z T 



C 1X T C 2 y 7 



z T W 2 z = 



C lX T C 2 y 7 



Al+a 2 B 12 Bl 2 a(A x B l2 + B l2 A 2 ) 
a(A 1 B 12 + B 12 A 2 f 



1 A 2 + a 2 Bf 2 B 12 



C x x 
C 2 y 



= C\x T A\x + C 2 2 y T A\y + a 2 (C\x T B X2 B\ 2 x + C 2 2 y T Bl 2 B 12 y) + 2aC 1 C 2 x T (A ± B 12 + B 12 A 2 ) y 
= C\X\ (Ax) + C 2 2 X\ (A 2 ) + 2C X C 2 6 + a 2 (c 2 \\B^x\\l + C 2 \\B 12 y\\fj 

where 9 — a (Ai (Ai) + Ai (A 2 )) x T B\ 2 y. Let Xl be the best possible lower bound z T W 2 z via the optimal linear 
combination (C\ and C 2 ) of x and y. Thus, 

X L = max C 2 X 2 (A x ) + C 2 X\ (A 2 ) + 2CxC 2 B + a 2 (c\ llS^Ha + C l H^llf) 
cf+c 2 =i \ ' 

We use the Lagrange multipliers method and define the Lagrange function as 

A = C 2 X 2 (A x ) + C 2 X 2 (A 2 ) + 2CxC 2 e + a 2 (c 2 \\Bj 2 xf 2 + C 2 \\B 12 y\\ 2 2 ) - » (C 2 + C 2 - l) 
where fi is the Lagrange multiplier. The maximum is achieved at the solutions of 

= 2Cx\ 2 (Ax) + 2aC 2 (Ai (Ax) + Xx (A 2 )) x T B 12 y + 2a 2 C\ \\Bf 2 x\\l - 2Cxv = 

oCx 
BA 

— = 2C 2 \\ (A 2 ) + 2aCx (Xx (Ax) + Xx (A 2 j) x T B 12 y + 2a 2 C 2 \\B 12 y\\l - 2C 2i i = 
uC 2 
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C\ + C| - 1 = 



which lead to ^CiJ^4- + C 2 -^^j /2 = Xl — M = 0. Hence, the maximum Xl is achieved at the solution of 

CxXj (Ax) + C 2 6 + a 2 Cx \\Bj 2 xf 2 - CxX L = 
C 2 X\ (A 2 ) + CxO + a 2 C 2 \\B 12 y\\ 2 2 - C 2 X L = 



that is 



This leads to 



det 



Xl (Ax)+a 2 \\Bl 2 x\\l-X L 



X 2 (A 2 ) + a 2 \\B 12 y\\ 2 2 - X L 



A 2 - (a 2 (Ax) + a 2 \\B^x\\l + X\ (A 2 ) + a 2 \\B 12 y\\f) A l + (a 2 (Ax) + 



= 



a 2 \\Bf 2 x\ h 



) (x\(A 2 ) + a 2 \\Bx 2 y\\ 2 2 )-e 2 = Q 
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Hence, 



(a? (A 1 ) + a 2 \\B? 2 x\\ 2 2 + X 2 (A 2 ) + a 2 \\B l2 y\\ 



X 2 (A,) + a? \\B? 2 x\\ 2 2 + A? (A 2 ) + a 2 \\B 12 y\\ 2 2 ) - 4 ((a? (A 1 ) + a 2 \\B^x\Q (a? (A 2 ) + « 2 U^H*) - 



X 2 (A 1 ) + a 2 WB&W; + X 2 (A 2 ) + a 2 \\B 12 y\\ 



which is obtained when 



\ 



'A? (A x ) + a 2 \\B{ 2 x\\: - X 2 (A 2 ) - a 2 \\B 12 y\ 



Ci = 



C 2 = 



A. 4 Proof of Theorem [9] 



9 2 + [X L -X 2 (Ai)-a 2 \\BT 



12*112 



X L -X 2 {A^-a 2 \\B( 2 x 



^je 2 + (x L -x 2 (Ai) 



a 2 \\ B i2 x \\ 2 



2N 



Any vector z of size 2N with zz T = m can expressed as a linear combination of the eigenvectors (zi, z 2 , z 2 m) 
of matrix W 

z 

where J2 2 =i c l = 1- Hence, 



= E C * 



z T W s z 



'2N \ /2N 

' CiZi I I Vr, H ;, 
=1 / \i=l 
T 



\i=l 

(2AT 
i=l 

2 N / 2 N . s 



E c ^ s 

Vi=l / 
/ 2JV \ 

E 04 



i=l Al . 



Hence. 



According to Lyapunov's inequality, 



'z T W k z s 
limk^oo [ — ^ — ) = Ai (W) 



(E[\X\ S ]) 1/S < [E \X\ l 



l/t 



when < s < t. Taking Pr[Af = j^-] = cf , we have 



2iV 



2iV 



i=l 1 i=l 



(E[\X\ S ]) 1/S <{E [\X\ 



l/k 



2N XI 



2^ 

*A* 



since fc is even and k > s > 0. 
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A. 5 Proof of Theorem [TD] 



Ai (W) — max 

x T x+y T y— 1 



max 

x T x-\-y T y— 1 



T T 

a; y 



(A + aB) 



x 
V 



T T 

x y 



A 



x 
V 



T T 

x y 



D 



x 

y 



< max (a; Aix + y A2y) + a max 

x T x+y T y— 1 x T x+y T y—l 



T T 
X y 



B 



x 
V 



= max(Ai(Ai),A 1 (A 2 )) +aX 1 (B) 
The inequality is due to the fact that the two terms are maximized independently. The second term 



Ai (B) = max (x T B 12 y + y T B^ 



is equivalent to the system of equations 



x T x+y T y— 1 

= 2 max x T B\ 2 y 

x T x+y T y— 1 



B 12 y = \ 1 {B)x 
B l2 y = X l {B)x 
x x + y y = 1 



or 



Bf 2 B 12 y = X^Bfy 
B 12 Bf 2 x = Ai(B) 2 .t 
x T x + y T y = 1 

which is to find the maximum eigenvalue (or more precisely the positive square root) of the symmetric positive 
matrix B\ 2 B\ 2 

Ai(£>) = , msix.x T B\ 2 Bi 2 x 

V x 2 — l 



This actually proves Lemma El the property Xi(B) = y/ X\(Bi 2 Bf 2 ) of a bipartite graph B. 

A. 6 Proof of Theorem [TTI 

The explicit expression up to the second order reads 

(A + aB) + azW + a 2 z (2) + Q{a ^ = ( A (o) + qA (i) + a 2 X (2) + ^ ^(o) + az (i) + Q 2 Z (2) + 0(q3) ^ 

(10) 

The zero order expansion is simply 

Az (0) =A (0) Z (0) 

The problem at zero order becomes to find the maximum of 

z (0)T Az (0) _ ( z (0) ; Az(°)) 



In the non-degenerate case, 



Hence, 



Z (0)T Z (0) ( Z (0) ;Z (0)) 

(zV^AzW) _ (x,A lX ) 
(z(°),z(°)) ~ (x,x) 

X^^X 1 (A 1 ) 
(z {Q) ) T = [x T ,0 T ] 



MM) 



1G 



where the first N elements of z^ are x and the rest N elements are all zeros. Let us look at the first order 
correction. Imposing the identity for the first order expansion in (|10[) gives 

Az^+Bz^ =X^z^+X^z^ (11) 

Furthermore, we impose the normalization condition to z (see ©), which leads to 

=0 

The first order correction to the principal eigenvector is orthogonal to the zero order. Plugging this result in 

CD 

(f>\ Az& + BzW) = A<°> (*<°>, + A (1) (>\ z {0) 
(A T z^\z^) + (z^\Bz^) =AW 



that is 



Since (z(°)) T = ( x T T ) and B 



U°\BzW)=\W (12) 
the first order correction is this non-degenerate case is 



B 12 
Bj 2 

null A^ = 0. Equation (fTTj) allows us to calculate also the first order correction to the eigenvector 



AzW+BzW=\(°hM 
(a - AW/) zW = -Bz<® 

(A — X^l) is invertible out of its kernel {A — X^ i) z = (that is the linear space generated by z^) and since 
Bz (o) ± z (o) we have 

«« = (AW/-A)~ 1 fl«W (13) 

Let us look for the second order correction. Imposing the identification of the second order term of (|T0l) we 
obtain 

Az^ + BzM = A(°M 2 > + X^z^ + A< 2 V°) 
Projecting this vectorial equation on z^ provides the second order correction to A 

(z(°\ Az^ + Bz W) = A<°>(z<°\ z (2) ) + A« (z<°>, z< x >) + X^ (*<°>, *<°) 
XW = ( z W,AzW) + ( Z W,BzW) - A(°)( z (°),z( 2 )) 
= A(°)(/°), z<») + ( z W,BzW) - A(°)(z(°), z^) 
= (zW,BzW) 



Substituting (fT5]l gives 



A (2) 



-,(0) 



(14) 



which can be further expressed as a function of the largest eigenvalue/eigenvector of individual network Ax, A 2 
or their interconnections B\2- Since 



Bz<® 



u Bxi 

>12 " 



bTo o 



X 
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wc have 



= ( (x^) T B 12 
(xM) T B 12 



-l 



(X^I-A^) WO \ 

v (A(°)/-A 2 ) J ^ Bf^W J 

' (X^I-A^ 1 \ / \ 

v (A(°)/-A 2 )- 1 J ^ ) 

which finishes the proof. 

A.7 Proof of Theorem [H 

In this case, the solution of the zero order expansion equation 

Az (0) =A (0) 2 (0) 

can be any combination of the largest eigenvector of the two individual networks x and y: 

= c\x + c 2 y 
cj+c 2 2 = l 

and X 1 - ' 1 — Ai {Ai) — Ai (A 2 ). The first order correction of the largest eigenvalue correction in the non-degenerate 
case (p~2]) holds as well for the generate case 



(z^\Bz^) =A« 



which is however non-zero in the degenerate case due to the structure of z^ and is maximized by the right 
choice of c\ and c 2 . Thus, 

X 1 (W)=max(x 1 (A 1 )+a(z(°\Bz(°y + o(a 2 ) 

= Xi(Ai) +maxacic 2 ((Bi 2 y,x) + (BT 2 x,y)) + o(a 2 ) 

ci,c 2 

= Ai(Ai) + -a ((B 12 y, x) + (Bj 2 x, y)) + o(a 2 ) 

= Ai(Ai) + -a (x, B 12 y) + o{a 2 ) 

where C\C 2 is maximum when c\ = c 2 = \/y/2. 

One may also evaluate the second order correction to the largest eigenvalues of the degenerate case. The 
following results we derived in the non-degenerate case hold as well for the degenerate case 

f A< 2 ) = (z(°\BzW) 

\ AzW+BzW = \<f>) z (i) +x {i) z {0) 

The latter equation allows to calculate the first order correction to the dominate eigenvector z^: 

i\Wl-A)zW = (B-\W)zW> (15) 

Any linear equation admits solutions when the constant term (i? — A^ 1 ') z^ ' is orthogonal to the kernel of the 
adjoint matrix of X^I — A. 

Ker{X {0) I-A) = \v: (X^I-A)v = o\ 
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Apart from pathological cases, each the two interactive nets have non degenerate maximum eigenvalues (A^ -*) 
corresponding to the dominant eigenvectors x^ and y^: 

\ A 2 yW = A(°y°); 

in this case, the kernel of the matrix X^I — A is just the linear space generated by the maximum eigenvalue of 
the single nets: 

{ by<"> J ' 

As stated, the constant term is orthogonal to the entire kernel: 



v(B- At 1 )) «(°) = ( ax® by<® ) ■ (b - A«) z<® = a {[B^x, y) - A«) + b ((*, B 12 y) - A«) 



0. 



Therefore the solution of the previous equation exists and all solutions differ by a vector in Ker(X^°'I — A). It 
is worth stressing that the value of A^ 2 -* does not depend on this extra terms, that is A*- 2 -* is invariant under the 
transformation: 

J xP) -> z« + ax<® 
\ V {1) -> V { 

providing the normalization condition is respected: 



,(!) -V „(1) -U6y(0); 



(z^fz^ = 

that is 

(x^) T x^ + (yWfyW = -> (z(°>) T (z« + a,' ') + (2/ (0) ) T (y (1) + &2/ (0) ) = 
that leads to a = b. Let us apply the transformation to A^ 2 ) : 

A(2) = J_ f fl„y(0) ^ f ^ ^ . (2) = J_ f Bla „(0) ^ f +oxW 



A< 2 ) = A< 2 > + a UyWfBTtxW - {x^f B 12 y^\ = A< 2 > 

Therefore we are allowed to select a definite solution as we where fixing a gauge. We will impose the orthogonality 
of x^ with x^ and j/ 1 ) with 

The linear operator X^I — A\ + x^ (x^) T ) and X^I — A\ operate identically over all vectors orthogonal 
to x^ (as all constant terms are in our case), while it behaves as an identity in the linear space generated 
by x^ \ Therefore to fix the gauge one may substitute X^I — Ay with X^I — A\ + x^ (x^) T ). The same 
argument holds for A 2 . This allows us to provide A^ 2 ^ an explicit expression and to calculate it algebraically. 
Equation [15] in components reads: 

J (AW/- AJx™ = B 12 yW - \W X (0) 
\ (\Wl-A 2 )yW = B&xW -AWyW; 

that, after fixing the gauge, becomes: 

J (A(°)J- A 1 +x^(x ( -°Y)x^ = fl 12 y(°) - AWorW 
\ (X^I-A 2 +y (a \y^) T )y {1} = Sf 2 ^°) - A^J/W; 

and hence the first-order correction to the dominant eigenvector can be 

J x<U = (X^I-A 1 +x(°\x^) T )- 1 (B 12 y(°)-X^x^) 
1 y W = (X^I - A 2 +y(°\y^) T r 1 (Bl 2 x^ - X^y^); 
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The second order correction can be finally be calculate algebraically, resorting to the spectral properties of 
the isolated networks, 

B ^y {0) V ( (A(°)/-A 1 + x(°)( a ;( )) T )- 1 o \( b 12 j/( 0) -a (1 ^ (0) \ 

~2\B{ 2 x^ J \ (X^I-A 2 +y^(y^) T r 1 ) \ Bj 2 x^ - \Wy<-°) J 

that is: 

A (2) = (y(0yf B T 2{x (0) I _ Ai+x (0)^(0) ) Tyl iBi2y (0)_ 

(16) 
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